module common_functions

	! self-defined functions based on matlab intrinsic functions
	implicit none
	real(kind=8), parameter :: pi = 3.141592653589793d+0
	real(kind=8), parameter :: cf_eps = 1.d-15
	integer, parameter :: no_digits = 20	! number of digits

	interface prtVec
		module procedure prtVecI, prtVecR, prtVecZ
	end interface prtVec

	interface prtMtx
		module procedure prtMtxI, prtMtxR, prtMtxZ
	end interface prtMtx

contains

    subroutine cf_bsxfunvm_d(op,v,m,rc,mm)
        ! op: 't' for times, 'd' for divide, 'p' for plus, 'm' for minus
        ! v: vector of size vn
        ! m: matrix of rows==vn when rc=='c', cols==vn when rc=='r'
        ! rc: 'c' means v is a column vector, 'r' means v is a row vector
        ! outputs:
        !   mm = v op m
        character, intent(in) :: op
        double precision, intent(in) :: v(:)
        double precision, intent(in) :: m(:,:)
        character, intent(in) :: rc
        double precision, intent(inout), optional :: mm(:,:)

        integer :: vn, sm, ii

        vn = size(v)
        sm = size(m)

        if ( present(mm) ) then
            if ( rc.eq.'r' ) then
                if ( op.eq.'t' ) then
                    forall (ii=1:vn)
                        mm(:,ii) = v(ii)*m(:,ii)
                    end forall
                end if
                if ( op.eq.'p' ) then
                    forall (ii=1:vn)
                        mm(:,ii) = v(ii)+m(:,ii)
                    end forall
                end if
                if ( op.eq.'m' ) then
                    forall (ii=1:vn)
                        mm(:,ii) = v(ii)-m(:,ii)
                    end forall
                end if
                if ( op.eq.'d' ) then
                    forall (ii=1:vn)
                        mm(:,ii) = v(ii)/m(:,ii)
                    end forall
                end if
            end if
            if ( rc.eq.'c' ) then
                if ( op.eq.'t' ) then
                    forall (ii=1:vn)
                        mm(:,ii) = v*m(:,ii)
                    end forall
                end if
                if ( op.eq.'p' ) then
                    forall (ii=1:vn)
                        mm(:,ii) = v+m(:,ii)
                    end forall
                end if
                if ( op.eq.'m' ) then
                    forall (ii=1:vn)
                        mm(:,ii) = v-m(:,ii)
                    end forall
                end if
                if ( op.eq.'d' ) then
                    forall (ii=1:vn)
                        mm(:,ii) = v/m(:,ii)
                    end forall
                end if
            end if
            return
        end if
        if ( .not.present(mm) ) then
            if ( rc.eq.'r' ) then
                if ( op.eq.'t' ) then
                    forall (ii=1:vn)
                        m(:,ii) = v(ii)*m(:,ii)
                    end forall
                end if
                if ( op.eq.'p' ) then
                    forall (ii=1:vn)
                        m(:,ii) = v(ii)+m(:,ii)
                    end forall
                end if
                if ( op.eq.'m' ) then
                    forall (ii=1:vn)
                        m(:,ii) = v(ii)-m(:,ii)
                    end forall
                end if
                if ( op.eq.'d' ) then
                    forall (ii=1:vn)
                        m(:,ii) = v(ii)/m(:,ii)
                    end forall
                end if
            end if
            if ( rc.eq.'c' ) then
                if ( op.eq.'t' ) then
                    forall (ii=1:vn)
                        m(:,ii) = v*m(:,ii)
                    end forall
                end if
                if ( op.eq.'p' ) then
                    forall (ii=1:vn)
                        m(:,ii) = v+m(:,ii)
                    end forall
                end if
                if ( op.eq.'m' ) then
                    forall (ii=1:vn)
                        m(:,ii) = v-m(:,ii)
                    end forall
                end if
                if ( op.eq.'d' ) then
                    forall (ii=1:vn)
                        m(:,ii) = v/m(:,ii)
                    end forall
                end if
            end if
        end if


    end subroutine cf_bsxfunvm_d




    subroutine cf_dsort(a,m)
        ! a stupid sort code
        implicit none
        double precision, intent(inout) :: a(:)
        integer, intent(in) :: m

        integer :: ii, jj(1)
        double precision :: dtmp

        do ii = 1, m-1, 1
            jj = minloc(a(ii:m))
            jj(1) = ii + jj(1) - 1
            dtmp = a(ii)
            a(ii) = a(jj(1))
            a(jj(1)) = dtmp
        end do
    end subroutine cf_dsort

	subroutine prtVecI(v, fmt)
		integer, dimension(:), intent(in) :: v
		character(len=*), intent(in), optional :: fmt

		character(len=no_digits) :: ns

		if ( present(fmt) .and. (fmt.ne.'*') ) then
			write(ns,'(I10)') size(v)
			write(*,'('//trim(ns)//trim(fmt)//')') v
		else
			write(*,*) v
		end if
	end subroutine prtVecI

	subroutine prtVecR(v, fmt)
		real(kind=8), dimension(:), intent(in) :: v
		character(len=*), intent(in), optional :: fmt

		character(len=no_digits) :: ns
		integer :: ii, sv

		sv = size(v)

		if ( present(fmt) .and. (fmt.ne.'*') ) then
			do ii = 1, sv, 1
				write(*,'('//fmt//'$)') v(ii)
			end do
			write(*,*)
		else
			write(*,*) v
		end if
	end subroutine prtVecR

	subroutine prtVecZ(v, fmt, trans)
		complex(kind=8), dimension(:), intent(in) :: v
		character(len=*), intent(in), optional :: fmt
		character, intent(in), optional :: trans

		character(len=no_digits) :: ns, s1, s2
		integer :: ii, sv

		sv = size(v)

		if ( present(trans) .and. (trans.eq.'T'.or.trans.eq.'t') ) then
			if ( present(fmt) ) then
				do ii = 1, sv, 1
					write(s1,fmt) dreal(v(ii))
					if ( dimag(v(ii)).ge.0.d0 ) then
						write(s2,fmt) dimag(v(ii))
 						write(*,'("'//trim(adjustl(s1))//'+'//trim(adjustl(s2))//'i  ")')
 					else
 						write(s2,fmt) -dimag(v(ii))
 						write(*,'("'//trim(adjustl(s1))//'-'//trim(adjustl(s2))//'i  ")')
					end if
				end do
			else
				do ii = 1, sv, 1
					write(s1,*) dreal(v(ii))
					if ( dimag(v(ii)).ge.0.d0 ) then
						write(s2,*) dimag(v(ii))
 						write(*,'("'//trim(adjustl(s1))//'+'//trim(adjustl(s2))//'i  ")')
 					else
 						write(s2,*) -dimag(v(ii))
 						write(*,'("'//trim(adjustl(s1))//'-'//trim(adjustl(s2))//'i  ")')
					end if
				end do
			end if
			write(*,*)
			return
		end if

		if ( present(fmt) .and. (fmt.ne.'*') ) then
			do ii = 1, sv, 1
				write(s1,fmt) dreal(v(ii))
				if ( dimag(v(ii)).ge.0.d0 ) then
					write(s2,fmt) dimag(v(ii))
 					write(*,'("'//trim(adjustl(s1))//'+'//trim(adjustl(s2))//'i  "\)')
 				else
 					write(s2,fmt) -dimag(v(ii))
 					write(*,'("'//trim(adjustl(s1))//'-'//trim(adjustl(s2))//'i  "\)')
				end if
			end do
			write(*,*)
		else
			do ii = 1, sv, 1
				write(s1,*) dreal(v(ii))
				if ( dimag(v(ii)).ge.0.d0 ) then
					write(s2,*) dimag(v(ii))
 					write(*,'("'//trim(adjustl(s1))//'+'//trim(adjustl(s2))//'i  "\)')
 				else
 					write(s2,*) -dimag(v(ii))
 					write(*,'("'//trim(adjustl(s1))//'-'//trim(adjustl(s2))//'i  "\)')
				end if
			end do
			write(*,*)
		end if
	end subroutine prtVecZ

	subroutine prtMtxI(mtx, fmt, trans)
		integer, dimension(:,:), intent(in) :: mtx
		character(len=*), intent(in), optional :: fmt
		character(len=1), intent(in), optional :: trans

		integer :: i, mn(2)
		mn = shape(mtx)
		if ( present(trans) .and. trans .eq. 'T' ) then
			do i = 1, mn(2)
				call prtVecI( mtx(:,i), fmt )
			end do
		else
			if ( present(fmt) ) then
				do i = 1, mn(1)
					call prtVecI( mtx(i,:), fmt )
				end do
			else
				write(*,*) mtx
			end if
		end if
	end subroutine prtMtxI

	subroutine prtMtxR(mtx, fmt, trans)
        ! such as prtmtxr(A,'f6.1','n')
		real(kind=8), dimension(:,:), intent(in) :: mtx
		character(len=*), intent(in), optional :: fmt
		character(len=1), intent(in), optional :: trans

		integer :: i, mn(2)
		mn = shape(mtx)
		if ( present(trans) .and. trans .eq. 'T' ) then
			do i = 1, mn(2)
				call prtVecR( mtx(:,i), fmt )
			end do
		else
			if ( present(fmt) ) then
				do i = 1, mn(1)
					call prtVecR( mtx(i,:), fmt )
				end do
			else
				write(*,*) mtx
			end if
		end if
	end subroutine prtMtxR

	subroutine prtMtxZ(mtx, fmt, trans)
		complex(kind=8), dimension(:,:), intent(in) :: mtx
		character(len=*), intent(in), optional :: fmt
		character(len=1), intent(in), optional :: trans

		integer :: i, mn(2)
		mn = shape(mtx)
		if ( present(trans) .and. trans .eq. 'T' ) then
			do i = 1, mn(2)
				call prtVecZ( mtx(:,i), fmt )
			end do
		else
			if ( present(fmt) ) then
				do i = 1, mn(1)
					call prtVecZ( mtx(i,:), fmt )
				end do
			else
				write(*,*) mtx
			end if
		end if
	end subroutine prtMtxZ


	subroutine cf_lglnodes(n,x,w,p)
		!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		!
		! lglnodes.m
		!
		! Computes the Legendre-Gauss-Lobatto nodes, weights and the LGL Vandermonde
		! matrix. The LGL nodes are the zeros of (1-x^2)*P'_N(x). Useful for numerical
		! integration and spectral methods.
		!
		! Reference on LGL nodes and weights:
		!   C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, "Spectral Methods
		!   in Fluid Dynamics," Section 2.3. Springer-Verlag 1987
		!
		! Written by Greg von Winckel - 04/17/2004
		! Contact: gregvw@chtm.unm.edu
		!
		!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		implicit none
		integer, intent(in) :: n
		real(kind=8), intent(out) :: x(n+1)
		real(kind=8), intent(out) :: w(n+1)
		real(kind=8), intent(out), optional :: p(n+1,n+1)

		real(kind=8) :: pp(n+1,n+1), xold(n+1)
		integer :: k


        ! 		x = real((/0:n/),8)
        forall (k=1:n+1)
            x(k) = k-1
        end forall
		x = x/real(n,8)
		x = cos(pi*x)

		xold = 2.d0

		do while ( maxval(dabs(x-xold)).gt.cf_eps )
			xold = x
			pp(:,1) = 1.d0
			pp(:,2) = x
			do k = 2, n, 1
				pp(:,k+1) = ( (2.d0*k-1.d0)*x*pp(:,k)-(k-1.d0)*pp(:,k-1) )/real(k,8)
			end do
			x = xold - ( x*pp(:,n+1)-pp(:,n) ) / ( (n+1)*pp(:,n+1) )
		end do
		w = 2.d0/(real(n*(n+1),8)*pp(:,n+1)**2)
		x = -x
		if ( present(p) ) p = pp
	end subroutine cf_lglnodes

	subroutine cf_linspace(x_start,x_end, n, val)
		real(kind=8), intent(in) :: x_start, x_end
		integer, intent(in) :: n
		real(kind=8) :: val(n)

		integer :: ii
		real(kind=8) :: h

		if ( n.le.1 ) then
			write(*,*) 'n should be great than 1 in cf_linspace.'
			stop 1
		end if

		h = ( x_end-x_start ) / real(n-1,8)
		do ii = 1, n, 1
			val(ii) = x_start + h*(ii-1)
		end do
	end subroutine cf_linspace

	recursive subroutine cf_qsort_i(A, sizeA, cmp)
  		! quick sort a integer vector into ascending order
    	integer, intent(inout), dimension(:) :: A
    	integer :: sizeA
    	interface
    		function cmp(x,y)
    			integer, intent(in) :: x, y
    			integer :: cmp
   		 	end function cmp
    	end interface

    	integer :: iq

    	if( sizeA>1) then
      		call cf_partition_i(A, iq, cmp)
      		call cf_qsort_i(A(1:iq-1), iq-1, cmp)
      		call cf_qsort_i(A(iq:sizeA), sizeA-iq+1, cmp)
    		endif
  	end subroutine cf_qsort_i

  	subroutine cf_partition_i(A, marker, cmp)
  		interface
    		function cmp(x,y)
    			integer, intent(in) :: x, y
    			integer :: cmp
    		end function cmp
   		end interface
    	integer, intent(inout), dimension(:) :: A
    	integer, intent(out) :: marker

    	integer :: i, j
    	integer :: temp
    	integer :: x      ! pivot point

   		x = A(1)
    	i = 0
    	j = size(A) + 1
    	do
      		j = j-1
      		do
        		if( cmp(A(j),x).le.0 ) exit
        		j = j-1
      		end do
      		i = i+1
      		do
        		if( cmp(A(i),x).ge.0 ) exit
        		i = i+1
      		end do
      		if( i<j ) then
        		! exchange A(i) and A(j)
        		temp = A(i)
        		A(i) = A(j)
        		A(j) = temp
      		else
      			if( i==j ) then
        			marker = i + 1
        			return
      			else
        			marker = i
        			return
        		end if
      		end if
    	end do
  	end subroutine cf_partition_i


end module common_functions
